Localized bases for kernel spaces on the unit sphere * 

t 

E. Fuselier* T. Hangelbroek § F. J. Narcowich 1 J. D. Ward'! 

G. B Wright** 

March 4, 2013 



Abstract 

Approximation/interpolation from spaces of positive definite or con- 
ditionally positive definite kernels is an increasingly popular tool for 
the analysis and synthesis of scattered data, and is central to many 
meshless methods. For a set of N scattered sites, the standard ba- 
sis for such a space utilizes N globally supported kernels; computing 
with it is prohibitively expensive for large N. Easily computable, well- 
localized bases, with "small- footprint" basis elements - i.e., elements 
using only a small number of kernels - have been unavailable. Working 
on S 2 , with focus on the restricted surface spline kernels (e.g. the thin- 
plate splines restricted to the sphere), we construct easily computable, 
spatially well-localized, small- footprint, robust bases for the associated 
kernel spaces. Our theory predicts that each element of the local ba- 
sis is constructed by using a combination of only 0((log A) 2 ) kernels, 
which makes the construction computationally cheap. We prove that 
the new basis is L p stable and satisfies polynomial decay estimates 
that are stationary with respect to the density of the data sites, and 
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we present a quasi-interpolation scheme that provides optimal L p ap- 
proximation orders. Although our focus is on S 2 , much of the theory 
applies to other manifolds - § d , the rotation group, and so on. Finally, 
we construct algorithms to implement these schemes and use them to 
conduct numerical experiments, which validate our theory for interpo- 
lation problems on § 2 involving over one hundred fifty thousand data 
sites. 
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1 Introduction 

Approximation/interpolation with positive definite or conditionally positive 
definite kernels is an increasingly popular tool for analyzing and synthesiz- 
ing of scattered data and is central to many meshless methods. The main 
difficulty in using this tool is that well-localized bases with "small-footprint" 
elements - i.e., elements using only a small number of kernels - have been 
unavailable. With this in mind, we have two main goals for this paper. 

The first is the theoretical development of small-footprint bases that are 
well-localized spatially, for a variety of kernels. For important classes of 
kernels on § 2 , the theory itself predicts that a basis element requires only 
0(log(iV) 2 ) kernels, where N is the number of data sites. 

Previous numerical experiments on data sets, with N on the order of a 
thousand, used ad-hoc techniques to determine the number of kernels per 
basis element. The predictions of our theory, on the other hand, have been 
verified numerically on §> 2 for data sets with over a hundred thousand sites. 

Our second goal is to show how to easily and efficiently compute these 
small-footprint, well-localized, robust bases for spaces associated with re- 
stricted surface-spline kernels on the sphere § 2 . The kernels in question are 
spherical basis functions having the form 

k m (x, a) := (-l) m (l - x ■ a)™" 1 log(l - x ■ a), (1) 

for m = 2, 3, . . . (cf. [T71 Eqn. 3.3]). The kernel spaces are denoted S m (E) - 
these are finite dimensional spaces of functions obtained as linear combina- 
tions of k m , sampled at some (finite) set of nodes 5 C § 2 , plus a spherical 
polynomial p of degree m — 1, i.e. X^ eS a§& m (-, £) + £>(•). The coefficients 
involved satisfy the simple side conditions given in [5} 

The Lagrange functions which interpolate cardinal sequences: X^(C) = 
Sg^, (6 5, form a basis for S m (S). Recently, it has been shown in [TT] . 
for restricted surface splines, as well as many other kernels, that these func- 
tions decay extremely rapidly away from £. Thus, {xd^es forms a basis 
that is theoretically quite good (sufficient to demonstrate that the Lebesgue 
constant is uniformly bounded, among many other things). However, de- 
termining a Lagrange basis function generally requires solving a full linear 
system with at least N := #H unknowns, so working with this basis directly 
is computationally expensive. In this paper we consider an alternative basis: 
one that shares many of the nice properties of the Lagrange basis, yet its 
construction is computationally cheap. 

Here is what we would desire in an easily computed, robust basis {^}^es 
for S m CB). Each basis function should be highly localized with respect to the 
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mesh norm h := maxg g = dist(x, £) of 3. Moreover, each should have a nearly 
stationary construction. By this we mean that each basis element is of the 
form X^eTfO where the coefficients A^ >rj and the degree 

m — 1 polynomial are completely determined by fc m and a small subset 
of centers Specifically, we wish to satisfy the following requirements: 

i) #T(0 = c(iV) 

ii) \k( x )\ ^ <r( r /h)> r '■= dist(x,£). 

where the number of points influencing each basis function c(N) is constant 
or slowly growing with N, and the function cr(-) decays rapidly - at an 
exponential rate a(t) < Ce _1/ '*' or at least at a fast polynomial rate cr(t) < 
C(l + The B-spline basis, constructed from the family of truncated 

power functions (i.e., using (x—y)™ in place of k m (x, y)), is a model solution 
to the problem we consider. 

Main results. The solution we present is to consider a basis of "local 
Lagrange" functions, which are constructed below in Section [3j It has the 
following properties: 

• Numerical Stability. For any J > 2, one can construct a numerically 
stable basis with decay a(t) < C(l + \t\)~ J . 

• Small footprint. Each basis function is determined by a relatively 
small set of centers: c(iV) < M (log AT) , where the constant M is 
proportional to the square of the rate of decay J: M a J 2 . 

• L p stability. The basis is stable in L p : sequence norms ||c||^ of the 
coefficients are comparable to L p norms of the expansion ^£e" c €^? - 

• Near-best approximation. For sufficiently large J, the operator 
Qsf = X^es /(£)&£ provides near-best approximation. 

Preconditioners. Over the years practical implementation of kernel 
approximation has progressed despite the ill-conditioning of kernel bases. 
This has happened with the help of clever numerical techniques like multi- 
pole methods and other fast methods of evaluation [U HI [12] and often with 
the help of preconditioners [31 EJ EJ [23] . Many results already exist in the 
RBF literature concerning preconditioners and "better" bases. For a good 
list of references and further discussion, see py. Several of these papers use 
local Lagrange functions in their efforts to efficiently construct interpolants, 
but the number of points chosen to localize the Lagrange functions are ad 
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hoc and seem to be based on experimental evidence. For example, Faul and 
Powell, in [7J, devise an algorithm which converges to a given RBF inter- 
polant that is based on local Lagrange interpolants using about thirty nearby 
centers. Beatson-Cherrie-Mouat, in [2], use fifty local centers (p. 260, Ta- 
ble 1) in their construction along with a few "far away" points to control the 
growth of the local interpolant at a distance from the center. In other work, 
Ling and Kansa p3| and co-workers have studied approximate cardinal basis 
functions based on solving least squares problems. 

An offshoot of our results is a strategy for selecting centers for precon- 
ditioning (as in [7J and [2]) that scales correctly with the total number of 
centers N. We demonstrate the power of this approach in Section [JJ where 
the local basis is used to successfully precondition kernel interpolation prob- 
lems varying in size by several orders of magnitude. 

Organization. We now sketch the outline of the remainder of the ar- 
ticle. In Section [2] we give some necessary background: in Section 2.1 



we 



treat analysis on spheres and in Section 2.2 we treat conditionally positive 
definite kernels. Section [3] presents the construction of the local Lagrange 
basis. Much of the remainder of the article is devoted to proving that this 
basis has the desired properties mentioned above. However, doing this will 
first require a thorough understanding of the (full) Lagrange basis {x^lfeSj 
which we study in detail in Sections [4] and [5j 

In Section [4] we consider the full Lagrange basis: the stable, local bases 
constructed in [11] . We first numerically exhibit the exponential decay of 
these functions away from their associated center. A subsequent experiment 
shows that the coefficients in the expansion X£ = S^es ^£,C^n('' have 
similar rapid decay. These numerical observations confirm the theory in 
Section [5j where it is proven that the Lagrange coefficients indeed decay 
quickly and stationarily with respect to h as £ moves away from £. 

Section [6] treats the main arguments of the paper. This occurs roughly 
in three stages. 

1. The decay of the Lagrange function coefficients indicates that trun- 
cated Lagrange functions \£ = X^er(£) ^,(^m(', C) wru be- satisfac- 
tory. But simply truncating causes the function to fall outside of the 
space (moment conditions for the coefficients are no longer satisfied), 
so it is necessary to adjust the coefficients slightly. The cost of this 
readjustment is related to the smallest eigenvalue of a certain Gram 
matrix: a symmetric positive definite matrix that depends on the set 



T(£) - this is discussed in Section 6.1 



2. Section 



6.2 



estimates the minimal eigenvalue of the Gram matrix, 
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which is shown to be quite small compared to tail of the coefficients. 
Although the resulting truncated, adjusted Lagrange function decays 
at a fast polynomial rate and requires 0((logiV) 2 ) terms, it is still 
unsuitable because its construction requires the full expansion. 

3. The local basis of Section [3] has coefficients sufficiently close to those 
of X£ to guarantee that it too satisfies the above properties. This, the 



main theorem and its corollaries are given in Section 6.3 

In Section [7] we demonstrate the effectiveness of using the local basis to build 
preconditioners for large kernel interpolation problems. 

Generalization to other manifolds/kernels. Finally, we note that 
many of the results here can be demonstrated in far greater generality with 
minimal effort: in particular, most results hold for Sobolev kernels on man- 
ifolds (as considered in [10J) and for many kernels of polyharmonic and 
related type on two point homogeneous spaces (as considered in [llj). To 
simplify our exposition, we focus almost entirely on surface splines on § 2 . 
We will include remarks discussing the generalizations as we go along. 



2 Background 

2.1 The sphere 

We denote by § 2 the unit sphere in M 3 , and by \i we denote Lebesgue 
measure. The distance between two points, x and £, on the sphere is written 
dist(x,£) := arccos(x • £). The basic neighborhood is the spherical 'cap' 
B(ct,r) := {x £ S 2 : dist(x, a) < r}. The volume of a spherical cap is 
fi(B(a,r)) = 2tt(1 — cosr). 

Throughout this article, 3 is assumed to be a finite set of distinct nodes 
on § 2 and we denote the number of elements in 3 by #S. The mesh norm 
or fill distance, h := /i(3,§> 2 ) := max xgS 2 dist(x,H), measures the density of 

3 in S 2 . The separation radius is q := | min^^ dist(£, £), where £, ( G 3, 
and the mesh ratio is p= := h/q. 

The Laplace-Beltrami operator and spherical coordinates. Given 
a north pole on S 2 , we will use the longitude 6± 6 [0, 2ir) and the colatitude 
82 € [0, 7r] as coordinates. The Laplace-Beltrami operator is then given by 

1 d . J)_ , 1 d 2 
sim^y W 2 Sm( 2 j W 2 + shr 2 ^) dej ' 

For each I S N, the eigenvalues of the negative of the Laplace-Beltrami 
operator, —A, have the form U£ := £(1 + £); these have multiplicity It + 1. 
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For each fixed £, the eigenspace Tie has an orthonormal basis of It + 1 
eigenfunctions, {Xiiu=—ii the spherical harmonics of degree 1. The space 
of spherical harmonics of degree I < a is U a = ©^ <cr Hi and has dimension 
(a + l) 2 . These and are the basic objects of Fourier analysis on the sphere. 
In order to simplify notation, we often denote a generic basis for IT^ as 
((j)j)j =1 ...( CT+ i)2. We deviate from this only when a specific basis of spherical 
harmonics is required. 

Covariant derivatives. The second kind of operators are the covariant 
derivative operators. These play a secondary role in this article - they are 
a construction useful for defining smoothness spaces and for proving results 
about surface splines, but they play no role in the actual implementation 
of the algorithms. We present some useful overview here - a more detailed 
discussion, where the relevant concepts are developed for Riemannian man- 
ifolds (including § 2 ) is given in |1U[ Section 2]. A more complete discussion 
is not warranted here. 

We consider tensor valued operators V m , where each entry (V m )(j 1) _ jifc ) 
is a differential operator of order m - each index ij = 1 or 2, corresponding 
to the variables x\ = 9\ and x<i = 02- They transform tensorially and 
are the "invariant" partial derivatives of order m. For m = 1, V is the 
ordinary gradient, but when m = 2, V 2 is the invariant Hessian. In spherical 
coordinates, V 2 it is a rank two covariant tensor with four components: 

r) 2 f r)f f) 2 f 

<W)., = (VVk. = - cot( S2 )|(). 

For a discussion of covariant derivatives on a general compact, C°° manifold, 
we refer the reader to JTUJ (2.7)]. 

Of special importance is the fact that at each point x E S 2 there is 
a natural inner product (•, •) on the space of tensors. The inner product 
employs the inverse of the metric tensor. For the sphere this is a diagonal 
2x2 matrix with entries: and g 1,1 (x) = (sin 2 (02)) 1 and g 2,2 (x) = 1. The 
general form of the inner product for tensors is 

<v"\f,v m 5 >, = £ r/f,)), jr s ( I )) il ,.,., i > i '( I )...^(4 

This gives rise to a notion of the pointwise size of the mth derivative at x: 

|V m /(x)| := V(V"\f,V"\/%. 
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Smoothness spaces. This allows us to construct Sobolev spaces. For 
each m and each measurable subset ft C § 2 , the L2 Sobolev norm is 



«T(n) : = ( E / |V fe /(z)| 2 d//(z) 



1/2 



2.2 Conditionally positive definite kernels and interpolation 

Many of the useful computational properties of restricted surface splines 
stem from the fact that they are conditionally positive definite. We treat 
the topic of conditional positive definiteness for a general kernel k. 

Definition 2.1. A kernel k is conditionally positive definite with respect 
to a finite dimensional space IJ if, for any set of N distinct centers S, the 
matrixK^ := C))^ * s positive definite on the subspace of all vectors 

a € satisfying X]geH = f or P ^ n ■ 

Let ((/?j)j e N be a complete orthonormal basis of continuous functions. 
Consider a kernel 

k(x,y) ■= ^2k(j)(pj(x)(pj(y) (2) 

with coefficients fe G ^(N) so that all but finitely many coefficients k(j) are 
positive. Then k is conditionally positive definite with respect to the (finite 
dimensional) space 77 := span(</?j | j G J), where J = {j \ k(j) < 0}. 
Indeed, 

E E = E e qwCOwCO = E^^ii^Hms) > 

provided Yl^ a ^fj(.0 = for j ^ J (i.e., satisfying /c(j) < 0). 

Conditionally positive definite kernels are important for the following 
interpolation problem. Suppose S C § 2 is a set of nodes on the sphere, 
/ : S 2 — > R is some target function, and /|„ are the samples of / at the 
nodes in E. We look for a function that interpolates this data from the space 

S(k,~) :=S(k,E,n) : = |e>^0,O I E>^(0 = 0, Vp G 77 J +77, 

Provided S C S 2 is unisolvent with respect to 77 (meaning that = for 
all £ 6 S implies that p = for any p G 77), the unique interpolant from 
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S(k, H) can be written 

s(-) = E a ^(''£) + E 

where the expansion coefficients satisfy the (non-singular) linear system of 
equations: 




where K s = i,j = 1, . . . , N, and $ = * = l,---,^, 

j G i7. This interpolant plays a dual role as the minimizer of the semi-norm 
| • \k induced from the "native space" semi-inner product 



Namely, it is the interpolant to f having minimal semi- norm \u\k = \J (u, u)^. 



3 Constructing the local Lagrange basis 

The restricted surface splines k m (see ([!])) are conditionally positive definite 
with respect to the space of spherical harmonics of degree up to m — 1, 
i.e. II m _i. The finite dimensional spaces associated with these kernels are 
denoted as in the previous section: 

= {Eces^M-, I E 5es o^(0 = o, v<^ g n m _!} + n m _{5) 

The goal of this section is to provide an easily constructed, robust basis for 
S m (S). The fundamental idea behind building this basis is to associate with 
each £ G 3, a new basis function that interpolates over a relatively small set 
of nodes a function that is cardinal at £. 

Specifically, let T(£) be the n <C N nearest neighbors to the node £, 
including the node £; see Figure [T] for an illustration. Then the new basis 
function associated with £ is given by 



&(•) = E ^.cM^o+E**.^")' 

C6T(0 J=l 



(6) 
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where cftj are a basis for the spherical harmonics of degree < m — 1. The 
coefficients A^s and c^j are determined from the cardinal conditions 

X £ (C) = l 1 lf C = ^ and V At C MC) = 0. (7) 

These coefficients can be determined by solving the (small) linear system 

where y £ represents the cardinal data and the entries of the matrix follow 
from ([3]). We call X£ a local Lagrange function about £. 



^.■.■..v;T\:^:::::V^: 
fc::v.v.'- {■■•••©■.•*/.i-.- :••.•>% 



1 1 



Figure 1: Illustration of the centers that make up the local Lagrange basis. 
The solid gray and black spheres mark the set of N nodes making up E. 
The solid black sphere with a circle around it marks the node £ where a 
local Lagrange function xs. is to be computed. The the solid black spheres 
enclosed in the dashed circular line mark the set of n = M(logN) 2 centers 
T(£) used to compute %£. For each £ E S, a similar set T(£) is determined 
for computing 

The new basis for 5 m (H) will consist of the collection of all the local 
Lagrange functions for the nodes in S. It will be shown in Section 6.3 that 



Localized Bases for Kernel Spaces 



11 



by choosing the number of nearest neighbors to each £ as n = M(logiV) 2 
will give a basis with sufficient locality. The choice of M is related to the 
polynomial rate of decay of \£ away from its center and a priori estimates 
are given for M in Section [6.3| However, in practice it will be sufficient to 
choose M by tuning it appropriately to get the desired rate of decay. 

The exact details of the algorithm for constructing this basis then pro- 
ceed as follows: For each £ G H 

1. Find the n = M(logN) 2 nearest neighbors to £, T(£). 

2. Construct X£ according to the conditions ([7]), which amounts to solving 
the associated linear system ^ and storing the coefficients Ag, c^. 

We note each set T(£) can be determined in O(logiV) operations by using 
a KD-tree algorithm for sorting and searching through the nodes S. After 
the initial construction of the KD-tree, which requires 0(N(logN) 2 ), the 
construction of all the sets T(£) thus takes 0(N(logN) 2 ) operations. 



Before continuing, we note that our main results, given in Theorem 6.5 
and its corollaries, depend heavily on properties that this local Lagrange 
basis inherits from the full Lagrange basis {x^lges- Thus, much of what 
follows is spent on developing a working understanding of the full Lagrange 
basis and its connections to the local Lagrange basis. Even though the local 
Lagrange basis is the focus of our work, we will delay any further mention 



of {x^l^es until Section 6.3 



4 The full Lagrange basis: numerical observations 

In this section we numerically examine a full Lagrange basis function x% and 
its associated coefficients for the kernel k2(x, a) = (1 — x-a) log(l— x-a), the 
second order restricted surface spline (also known as the thin plate spline) 
on S 2 . First, we demonstrate numerically that x% decays exponentially away 
from its center. Secondly, we provide the initial evidence that the Lagrange 
coefficients decay at roughly the same rate, which is proved later in Theorem 
PI 

The full Lagrange function centered at £ takes the form 

where is a degree 1 spherical harmonic. In this example, we use the 
"minimal energy points" of Womersley for the sphere - these are described 
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Geodesic Distance From Center 

Figure 2: Maximum latitudinal values of the Lagrange function for the ker- 
nel k2(x, a). This experiment was carried out in double precision arithmetic, 
and the plateau at roughly 10 -11 occurs due to ill conditioning of the collo- 
cation matrices and truncation error. 



and distributed at the website [25] jj Because of the quasi- uniformity of the 
minimal energy point sets, it is sufficient to consider the Lagrange function 
centered at the north pole £ = (0, 0, 1). 

Figure [2] displays the maximal colatitudinal value^] of \xd\- Until a 
terminal value of roughly 10~ n , we clearly observe the exponential decay of 
the Lagrange function, which follows 

Mx )\<C L e^(-u L ^^\ (9) 

(this "plateau" at 10 -11 is caused by roundoff error - see Figure [4]). The 

1 These point sets are used as benchmarks: each set of centers has a nearly identical 
mesh ratio, and the important geometric properties (e.g., fill distance and separation 
distance) are explicitly documented. 

2 The function \i ls evaluated on a set of points (01,62) with no equispaced longitudes 
6\ 6 [0, 2?r] and m equispaced colatitudes 62 £ [0,7r]. 
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o N=900 
a N=2500 
□ N=10000 
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Geodesic Distance From Center 



Figure 3: Plot of coefficients for a Lagrange function in the kernel space 
S(k2, 2). This experiment was carried out in double precision arithmetic. 



estimate Q has in fact been proven in Theorem 5.3], where this and 
other analytic properties of bases for S m (E) were studied in detail. By fitting 
a line to the data in Figure [2] where the exponential decay is evident, one can 
estimate the constants vl and Cl, which in this case are quite reasonable. 
For example, the value of ul, which measures the rate of exponential decay, 
is observed to be close to 1.35 (see Table [I]). 

We can visualize the decay of the corresponding coefficients in the same 
way. We again take the Lagrange function centered at the north pole: for 
each £' 6 S, the coefficient l-A^'l in the expansion X( = Y1 -^€,C^('> + * s 
plotted with horizontal coordinate dist(£, £'). The results for sets of centers 
of size N = 900, 2500 and 10000 are given in Figure [3j The exponential 
decay seems to follow 



1-4 



<c, 



q exp 



dist(£,C) 



h 



Indeed, this is established later in Theorem 5.3 As before, we can estimate 



the constants v c and C c for the decay of the coefficients. Comparing Figures 
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[2] and [3j we note that the coefficient plot is shifted vertically. This is a 
consequence of the factor of q~ 2 in the estimate (15) below. Table gives 



estimates for the constants v c and C c , along with the constants involved in 
the decay of the Lagrange functions. 



N 


h x 


PX 




C L 




Co 


400 


0.1136 


1.2930 


1.1119 


0.8382 


1.0997 


0.5402 


900 


0.0874 


1.5302 


1.3556 


1.0982 


1.3445 


0.7554 


1600 


0.0656 


1.5333 


1.3513 


1.2170 


1.3216 


0.5946 


2500 


0.0522 


1.5278 


1.3345 


0.9618 


1.3117 


0.5494 


5041 


0.0365 


1.5304 


1.3395 


1.1080 


1.3158 


0.6188 


10000 


0.0260 


1.5421 


1.3645 


1.1934 


1.3369 


0.7291 



Table 1: Estimates of the decay constants v and C for Lagrange functions 
and coefficients on the sphere using the kernel k2(x, a), with relevant geo- 
metric measurements of the minimum energy node sets used. 



The perceived plateau present in the Lagrange function values as well 
as the coefficients shown in Figures [2] and [3] is due purely to round-off error 
related to the conditioning of kernel collocation and evaluation matrices. 
These results were produced using double-precision (approximately 16 dig- 
its) floating point arithmetic. To illustrate this point, we plot the decay rate 
of the Lagrange coefficients for the 900 and 1600 point node sets as com- 
puted using high-precision (40 digits) floating point arithmetic in Figure |4j 
The figure clearly shows that the exponential decay does not plateau, but 



continues as the theory predicts (see Theorem 5.3). 



5 Coefficients of the full Lagrange functions 

In this section we give theoretical results for the coefficients in the kernel 
expansion of Lagrange functions. In the first part we give a formula relating 
the size of coefficients to native space inner products of the Lagrange func- 



tions themselves (this is Proposition 5.1). We then obtain estimates for the 
restricted surface splines on S 2 , demonstrating the rapid, stationary decay 
of these coefficients. 
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Lagrange Coefficient Decay (40 digit arithmetic) 
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Figure 4: Size of the coefficients for a Lagrange function in the kernel space 
S(k2, 2). This experiment was carried out in Maple with 40 digit arithmetic. 



5.1 Interpolation with conditionally positive definite kernels 

In this section we demonstrate that the Lagrange function coefficients 
can be expressed as a certain kind of inner product of different Lagrange 
functions and Xc- Because this is a fundamental result, we work in 
generality in this subsection: the kernels we consider here are conditionally 



positive of the type considered in Section |2.2 

When u,v £ S(k, H) - meaning that they have the expansion u = 
X^es a uM'> O+Pu and v = X^ e5 a 2 ^k(-, £)+p v with coefficients (a^^gn _L 
(77)|= for j = 1,2 - then the semi-inner product is 



(«• v )k = ( Y °u fc (-» o> Y a 2,^(-' o ) = Y Y a u a 2,cH^ 0- 

(This follows directly from the definition Q coupled with the observation 
that for j <£ J, u{j) = E^h a i,^0')<Aj(6 and v(j) = J2^=. a uHj) ( t ) j(0-) 
We can use this expression of the inner product to investigate the kernel 
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expansion of the Lagrange function. 



Proposition 5.1. Let &(•,£) = ^2j £ ^k(j)(j)j(-)(j)j(^) be a conditionally pos- 
itive definite kernel with respect to the space LI = span„ e j- <pj, and let H be 
unisolvent for LI . Then Xri £ S(k, S) (the Lagrange function centered at rj) 
has the kernel expansion Xvi x ) = X^es ^ri,^(x, £) + P( with coefficients 

A ri = (Aj^)eeB = «X<(z)>X»/(z)>fc) SgS - 

Proof. Select two centers (, rj € S with corresponding Lagrange functions xc 
and x»7 £ 3). Denote the collocation and auxiliary matrices, introduced 
in Section 2.2 by K= = C))^ an d $ = j- Because and 

are both orthogonal to (77) |~, we have 

(X{,Xy)k = Yl A C,ti A v,b k (ti,&) = (KhA^A^^h). 

Now define P := $($*$)- 1 $* : £ 2 (E) -> (IIj-JIs C 4(5) to be the 
orthogonal projection onto the subspace of samples of 77 on S and let 
P -1 - = Id — P be its complement. Then for any data y, ^ yields co- 
efficient vectors A and c satisfying P A = A and P ± <frc = 0, hence 
P- L K a J J - L A = P-LRhA = P ± y. Because P x : £ 2 (a) -> 4(H) is also an 
orthogonal projector, and therefore self-adjoint, it follows that 

(xd x )'Xr,i x ))k = (K H A C , A,,)^) 

= (P^KbA^A^.^b) 
= {P x e A v ) MS) . 

In the last line, we have introduced the sequence = (<$c,f)£ss for which 
K=A^ + P(\e = e ( which implies that P-'-KbA^ = P^e^. Using once more 
the fact that P ± is self-adjoint, and that A^ is in its range, we have 

{xd x ),Xy( x )}k,j = {P ± e A v ) = {e ( ,P ± A ri ) = (e c ,A„) 

and the lemma follows. □ 

The next result involves estimating the norms ||a||^ 2 (H) an d ll c ll£ 2 (J')' 
where a and c are as in ([3]). It will be useful later, when we discuss local 
Lagrange functions. The notation is the same as that used in the proof 
above. In addition, because A; is a conditionally positive definite kernel for 
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II, the matrix P K=P is positive definite on the orthogonal complement 
of the range of <£. We will let $ be the minimum eigenvalue of this matrix; 
that is, 

<?:= min (P ± K s P ± a, a) > 0. 

||P x a||=l 

Proposition 5.2. Suppose a and c satisfy Lei G= = <&*<!>. T/ien, 
\\ah 2 (E) < ^'llylli^B) < ^ i y#H||y||^ (S ) 

and 

\\c\\ MJ) < 2||A;|| 00 ||G s 1 || 1 /2 l9 -i #H || y || £2(H) 

< 2||A;|| 00 ||G s 1 || 1 / 2 1 9- 1 (#H) 3 / 2 ||y||, oo(H) . 



Proof. From (J3J) and the fact that P 1 - projects onto the orthogonal comple- 
ment of the range of <3?, we have that P J "K=P- L a = P^y and that P ± & = a. 
Consequently, 

tf||a||f 2(3) = ^||P ± a||| 2(3) < (P ± K s P ± z, a ) < ||a||, 2(H )||P ± y||^( H ). 

The bound on ||a||^ 2 (o) follows immediately from this and the estimate 
llylk(H) < V^^llyllfooCH)- To get the bound on ||c||^( 3); note that <£>c = 
Py — PK=a and, hence, that 

ll*c||< a(a) < ||Py[| /a(s) + UPKaPJ-all^s) 

< IIPyllfe^+^lPKHP^IHIPVlk^) 

We also have that ll^cllf^s) = (<3?*<I>c,c) > \ min (®*&)\\c\\ 2 e ^jy However, 
A m m(^*$) = IK^*^) -1 !! -1 ! which implies that 

\\c\\ HJ) < IK***)- 1 !! 17 '^^® = Wa^n**^ 

Next, note that the following inequalities hold: UPKgP^H < ||Kh|| < 

#311*1100, ||Py||, 2 ( H ), ll^yll^cs) < l|ylk(3) < VWhW^E), and > 

1. Applying these to the inequality 

\\4hj) < ||G H 1 || 1 / 2 (||Py||, 2(s) + #3||A;|| 00 ^ 1 ||PVll£ 2 (H)) 
then yields the desired bound on ||c||£ 2 (j), completing the proof. □ 
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5.2 Estimating Lagrange function coefficients 

In [TT], it has been shown that Lagrange functions for restricted surface 
splines decay exponentially fast away from the center. We can use these 



decay estimates in conjunction with Proposition 5.1 to estimate the decay 
of the coefficients 

Recall that the eigenvalues of —A are Xe = £(£+1). Let Q(z) := U™ =1 (z— 
\ u -i) = J2T=i ^vZ v ■ The kernel k m : E 2 x S 2 — >• E has the expansion 

oo I 

k m (x,a) =J2%{i) Y l( x ) Y l(<x)i 

where, for £ > m, k{£) = C m Q(X e )-\ with C m = 2 m+1 7rT(m) 2 H3 Eqn. 3.3]. 
From the expansion, one sees that k m is conditionally positive definite with 
respect to II m _i. Kernels such as k m are said to be of polyharmonic or related 
type; they have been studied in [11]. The kernel k m acts as the Green's 
function for the elliptic operator C m := C m l Q{— A) (cf. [11] Example 3.3]), 
in the sense that 

/=/ k m (-,a)£ m [f(a) -p f (a)]da+p f , 
Jm 

where pf is the orthogonal projection of / onto Tl m -i- 

The native space "inner product" on subsets. In [11] it was shown 
that for any k G N, the operator (V fc )*V fc (which involves (V fc )* the adjoint - 
with respect to the L2(S 2 ) inner produc t - o f the covariant derivative opera- 
tor V k which was introduced in Section 2.1 ) can be expressed as Y2v=o d v ky 



with dk = (— 1) . Consequently, any operator of the form Y^j=o Cj(V- 5 )*V- 7 
can be expressed as Y^=odv^ u with dk = (— l) fc Cfc and vice-versa: 

V(rfo, . . . d m ) 3(c , . . . , c m ) with d m = (-l) m c m 

m m 

and J^A" = ^ c i( vi )* vj - ( 10 ) 
u=0 j=l 

Because C m = (7~ 1 Q(-A), it follows that C m = Ejlo c i( vi )* vi > with 
c m = C™ 1 , an d so the native space semi-inner product, introduced in ([I]), 
can be expressed as 

{u,v) km = {C m u,v) L2i p ) = \ (3(u,v) x du.(x) 

Js 2 
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with (3(u,v) x = YlT=o c k(^ u,X7 v) x and co,...,c m are the appropriate 



constants guaranteed by (5.2). The latter expression allows us to extend 
naturally the native space inner product to measurable subsets Q of S 2 . 
Namely, 



(u,v)n,k m ■= / /3{u,v) x dn.{x). 
Jn 

This has the desirable property of set additivity: for sets A and B with 
fi(AnB) = 0, we have {u, v) A uB,k m = {u, v) A ,k m + {u, v) B ,k m - Unfortunately, 
since some of the coefficients c& may be negative, f3(u,u) and {u,u)Q^ m 
may assume negative values for some u: in other words, the bilinear form 
(u, v) h- > (u, w)n,fc m is only an indefinite inner product. 

A Cauchy-Schwarz type inequality. When restricted to the cone of 
functions in W^"(f2) having a sufficiently dense set of zeros, the quadratic 
form {u,u)fi t k m is positive definite. We now briefly discuss this. 

When Q has Lipschitz boundary and u has many zeros, we can relate the 
quadratic form |||m|||q km := (u, u)Q^ m to a Sobolev norm IMI^m^- Arguing 
as in [TTJ (4.2)], we see that 



< / fi(u,u) x du.{x) 
Jn 

< {f^\ c j\)\\ u \\w^(ny 



If u\~, = on a set 3 with /i(H, SI) < ho with ho determined only by the 
boundary of Q (specifically the radius and aperture of an interior cone con- 
dition satisfied by <9S7), Theorem A. 11 of [llj guarantees that ||u|| 2 ^ /m _i^^ < 

Ch 2 \ 

u \w™(Q) w ^h C depending only on the order m and the roughness of the 
boundary (in this case, depending only on the aperture of the interior cone 
condition). Thus, by choosing h < h* , where h* satisfies the two conditions 

h* < h and C{h*f x (max|c,|) < (11) 

j<m 2 



we have 



i (f2 ) < IHIn,fc m ^ (™^' C J'I ) H"liir.;"(<>r 



The threshold value h* depends on the coefficients Cj as well as the radius 
Rci and aperture of the cone condition for f2. When Q is an annulus of 
sufficiently small inner radius, the cone parameters can be replaced by a 
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single global constant, and h* can be taken to depend only on cq, . . . , c m . In 
other words, only on k m - cf. [Ill Corollary A. 16]. 

A direct consequence of this is positive definiteness for such functions, 
I u I q km > with equality only if u\n = 0. From this, we have a version 
of the Cauchy-Schwarz inequality: if u and v share a set of zeros Z (i.e., 
u\z = v\z = {0}) that is sufficiently dense in Q, then 

\(u,v)n )km \ < |H|| Q)fcm Hn ifcm (12) 

follows (sufficient density means that h(Z,£l) < h* as above). 

Decay of Lagrange functions. |lll Lemma 5.1] guarantees that the 
Lagrange function x% satisfies the bulk chasing estimate: there is a fixed 
constant < e < 1 so that for radii r the estimate 

IIXellwy»(B«K,r)) < 4xdw^{B-^,r-^)) 

holds. In other words, a fraction (roughly 1 — e) of the bulk of the tail 
\\X£ \\w^-{B c {i,r)) is to be found in the annulus £>(£, r) \ B(^, r — of width 
2^- oc h (with a constant of proportionality ^ that depends only on m). 
For r > 0, it is possible to iterate this n times, provided < r. It follows 
that there is v = — 4/io log e > so that 

llxdlw/ 2 m (B c (£,r)) < en ||xdlw 2 m (§ 2 ) - Ce~ ur ^ h \\x^\\w^(s 2 )- 
By [IH (5.1)] Qwe have 

\\Xt\\wp(B°(e,r)) < Cq l - m e-^. (13) 
This leads us to our main result. 

Theorem 5.3. Let p > be a fixed mesh ratio. There exist constants h* , 
v and C depending only on m and p so that if h < h* , then the Lagrange 
function xc, = Sges ^C?^™('> +P( ^ S m (E) has these properties: 

dist(a;, £) 
' h 
„2-2m_/' ,, dist (C,C) 



| XC |<Cexp(-^ (14) 



|4^| < Cq 2 ' 2m exp [^-v ^ J ■ ( 15 ) 
ci</ 2/ la|| £p(H) < || ^aadL p (§2) ^ c 2 q 2/p \\a\\ ep{s) . (16) 



3 This is simply a comparison of \£ to a smooth "bump" 0^ of radius q - also an 
interpolant to the delta data (8^), but worse in the sense that Ixdllfc — lll'/'lllfe ~ this idea 
is repeated in the proof of Theorem 



5.3 
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Proof. The bounds ( 14 ) and ( 16 ) are given in Theorems 5.3 & 5.7] . Only 



(15) requires proof. By Proposition |5 . 1 1 and set additivity, we have that 

where we employ the hemispheres: fi( = {a 6 S 2 | dist(a, £) < dist(a, £)} , 
Og = {a G § 2 | dist(a, £) < dist(a, £)} • Modulo a set of measure zero, Q,£ = 

s 2 \n ( . 



We apply the Cauchy-Schwarz type inequality (12) to obtain 



\ A C,t\ < \lxd\n okm \lxdn okm + \h(\in s ,k m M %fcm 



Since ft c C B c ((,r) := § 2 \ B ((, |dist(£,C)) and C B c (^r), we can 
again employ set additivity and positive definiteness (this time |||x^||q c k — 

III X£ III §2 k > wm ch follows from the fact that § 2 = f^Uf^ and that xs, vanishes 
to high order in - the same holds for X() to obtain 

\A U \ < Jmax\cj\ (llxcllw 2 -(i?-(C,r)) Ixdkm + ML WxnW wp(B^,r))) ■ 

The full energy of the Lagrange function can be bounded by comparing 
it to the energy of a bump function - for this is (j)^, which can be defined 
by using a smooth cutoff function a. In spherical coordinates (colatitude, 
longitude) around £, (f)^(8,(p) = a(9/q). This is done in [Til (5-1)] and we 
have that |||x£||| fc and Ixcllfc are bounded by Cq l ~ m . 



On the other hand, we can employ (13) to treat w^ l (B c ((.r)) and 
lx?lliy 2 m (s c (C,r)), which gives 

llx?llvK 2 m (B c (C,r))) llxellw 2 m (B c (C,r)) < Cq m e u h=Cq m e 2h . 



The bound (15) follows immediately from this. □ 



Remark 5.4. Because the proof doesn't really depend on S 2 , a nearly iden- 
tical proof works for any of the kernels with exponentially decaying Lagrange 



functions considered in [10(. II 1| . Specifically, we have this: Theorem 5.3 holds 
for compact, 2-point homogeneous spaces with polyharmonic kernels satisfy- 
ing £ m _L IT (cf. FIT]/ ) and for any compact, C°° Riemannian manifold, with 
the kernels being the Sobolev splines given in [10]. 
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6 Truncating the Lagrange basis 

We now discuss truncating the kernel expansion Lagrange function X£ = 
X<s ^£,C^"('' C) + Pi £ £m(2), replacing it with an expansion of the form 

X(= Yl h( k m(;C)+P^ S m (E), (17) 

C6T(0 

where T(£) C S is a set of centers contained in a ball B(£,r(h)) centered 
at £, where r(/i) and the A;,c' s wm be determined by A^, with £ £ T(£). 
We also assume that £ G T(£). Finally, to avoid notational clutter, we will 
simply use T rather than T(£). 



Our goal is to show that if x§ satisfies the properties (14), (15), and 



(16), then we may take r(h) = Kh\ log(/t)|, with K = K{m) > 0, while 



maintaining algebraic decay in h of the error \\x^ — Xclloo- F° r this choice 



of r(/t), a simple volume estimate (given at the end of Section 6.3) shows 
that the number of terms required for xs, is just C((log iV) 2 ) <C N, far fewer 
than the iV needed for xs,- 

Simply truncating at a fixed radius r(h) = Kh\ \og(h)\ is not suitable, 
however, because the truncated function X(, wm no longer be in the space 
5 m (S) (and thus {x§} wu l n °t act as a basis). To treat this, we must slightly 
realign coefficients to satisfy the moment conditions. 

A remark before proceeding with the analysis: Finding Xi m the way 
described below requires knowing the expansion for xe, an d carrying out the 
truncation above. This is expensive, although it does have utility in terms 
of speeding up evaluations for interpolation when the same set of centers is 
to be used repeatedly. The main point is that we now know roughly how 
many basis elements are required to obtain a good approximation to X(- 
The question of producing a good basis efficiently is left to the next section. 



6.1 Constraint conditions on the coefficients 

We would like X£ to be in the space S m (E), and so the -A^'s have to satisfy 
the constraints in the system Q: 

]T A^cMO = 0, j G J := (1, • • - , m 2 ), (18) 
Cex 

where is an orthonormal basis for Il m _i. Since the original x^'s are 

in S m (H), the A;,c' s m their expansions satisfy the constraint equations in 
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([3]). Splitting these equations into sums over T and its complement in 3 
and manipulating the result, we see that 

A=,cW(0 + °j> where a i Yl A uh(0, j e J. (19) 

The way that we will relate the two sets of coefficients is to define the vector 
(A^^gx to be the orthogonal projection of (A^^)^ e y onto the constraint 
space, which is the orthogonal complement of span{(/>j| T ,j < m 2 }, in the 
usual inner product for £2 (T). The equations below then follow: 

(Af,c)ceT - (^,c)cex = J^Tj^jW G s P an {(0i(C))ceT, j < ™ 2 } 
II A.C)C6T " (^,c)ceT|| £ 2 2{T ) = T-*G T r, [G T ] fcj := EceT&(0^(C), 

(20) 

where r is a column vector having the r^'s as entries. Let a be a column 
vector with the Oj's as entries. From the first equation above together with 
equations (18) and (19), r and a are related by a = G^t. If we make the 
rather mild assumption that T is unisolvent for the space Il m _i, then we 
can invert Gy. r = G^a, thereby obtaining t*G^t = a*GZ a. Using this 
in (20) and applying Schwarz's inequality, we obtain this bound: 



||E C6 r(^,C-A^)M-,C)|L < V# T H^Hz IIM00IMI2, (21) 

which we will make use of to establish the estimates below. 

Proposition 6.1. Assume that T is unisolvent for H m -i and that \\G^ W2 = 
0(| log h\~ 2 h~ 2 ^), for some fi > 0. If we take r{h) = Kh\\og{h)\, where K 
is chosen so that J := Kv — 2m — fi > then for h sufficiently small, 

\\xt-Xt\\°o<Ch J (22) 
\Xt(x)\ <C{1 + dist(x, 0/h)' J (23) 

Furthermore, when J > 2, the set {xs,} is L p stable: there are Cx,C 2 > 
for which 

Cig 2/p ||a||, p(H) < HEees^Xell^dp) < C 2 q 2 n*h v{ ~). (24) 
Proof. From (Il5]) and N < 4ir/vol(B(£, q)) < Cq~ 2 , we have that 



\ A t,c\ = 0(Nq 2 - 2m exp(-ur(h)/h)) < Ch Ku - 2m . (25) 
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Kv-2m 



Applying it to the a/s defined in (19) results in 1 1 c" 1 1 2 < Ch 
this in connection with {{Gy 1 ^ = 0{\ log /i|~ 2 hr 2 * 1 ), Q and 

% - X£ = $^(4e,< - A u) k m{-,0 ~ Ec^T^,C fc m(-, C), 

CeT 



Using 



yields (22). Next, from (14) we have 



\Xi\ < C ex P 



dist(x, £) 



< Cexp \-Kv 



< C 1 + 



dist(x,£) 



dist(x, £) 



Combining this with (22), using J = — 2m — \i > and manipulating, 
we arrive at (23). 

It remains to demonstrate the L p stability of (x^) for 1 < p < 00. When 
p = 1, we consider a sequence a = (ag)^ g = £ ^i(H). Let s := ^a^x$ an d 
s := E a £X£- From Holder's inequality and (22), we have ||s — < 
C|l a ll£i(H)^ J an d 

\\s ~ s|| ioo ( S 2) < C||a|| £oo(s) Efeslxs^) - xd x )\ ^ C|l a ll^c(H)/i V 2 - 

v -v ' 

<ATmax e ||x«-X£lLoo(S 2 ) 

Interpolating between these two inequalities - i.e., interpolating the finite 
rank operator a 1— >■ (s — s) - gives 



< C^-V^llall^p). (<f 



h~ 



After some manipulation, this bound and (16) imply that 



ciq 



2/PI 



p(3) (l - C^- 2 ) < ||5|| ip(s2) < c 2 ^||a||, p(H) (l + C/i 



2/PI 



Choosing /i so that Ch < 1/2 and letting Ci = ci/2 and C2 = 3c2/2, we 
obtain (pi). □ 



Remark 6.2. When there are no constraint conditions on the coefficients, 
this result holds for any of the strictly positive definite kernels mentioned 
in Remark 5.4, In particular it holds for Sobolev splines on a compact C°° 
Riemannian manifold. 
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6.1 



arc 



6.2 Norm of the inverse Gram matrix 

We now demonstrate that the conditions on G^ 1 in Proposition 
automatically satisfied. We will state and prove the results below for caps 
on S d , rather than just § 2 . Also, It is more convenient to use with II ^ rather 
than II m _i, because Tfi is notationally tied to the polyharmonic kernels k m 
as well as the spherical harmonics on § 2 . That said, we begin with the 
lemma below. 

Lemma 6.3. Suppose that S r := B(£, r) C S d is a cap of fixed radius r < ir, 
and that f C S r is finite and has mesh norm h<g := hs ry c g- In addition, let 
L > be a fixed integer and take Ul to be the space of all spherical harmonics 
of degree at most L. Then, there exists a constant cq := co(d,L) > such 
that when h<# < c^r we have 

W(0? > KSr)- 1 [ \^x)\ 2 dfi{x), for all <p G U L . (26) 

Moreover, the set 'if is unisolvent forU^- Finally, for every basis for IL^ the 
corresponding Gram matrices G<g and Gs r , relative to the inner products on 
l 2 (^€) and S r , respectively, satisfy 

\\G^h<KSr)\\Gs%. (27) 

Proof. Since since <p(x) and Tp{x) are spherical harmonics in IJ^, their prod- 
uct is a spherical harmonic of degree at most 2L. Thus, applying the 
nonnegative- weight quadrature formula in \16\ Theorem 2.1] to spherical 
harmonics of order 2L yields 



J> C |^(C)| 2 = f W(x)\ 2 d^{x), 



Since < < X^e 1 ^ w ( = M'SV), we have the inequality YlteV ^cl^OI 2 — 
^(Sr) Y^ceC I V^CO 1 2 - The inequality (26 ) follows immediately from the quadra- 



ture formula. To prove that C is unisolvent, suppose that ip 6 II ^ vanishes 



on C. By (26), we have that f s \ip(x)\ 2 du.(x) = 0. Since (p is in IL^, it is a 
polynomial in sines and cosines of the angles used in the standard param- 
eterization of S d , with ^ being the "north" pole. As a consequence, it is 
continuous on S r and, because J s \ip(x)\ 2 dfi(x) = 0, it is identically on 
S r . Finally, as a function of the angular variables in the complex plane, it 
is analytic, entire in fact, and can be expanded in a power series in these 
variables. The fact that it vanishes identically for real values of the angular 
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variables is enough to show that the coefficients in the series are all zero. 



Hence, ip = on S and is unisolvent for LT^. To establish (27), note 



that (26) implies that G<g — jJ>(S r ) 1 Gs r is positive semi definite. From the 



Courant-Fischer theorem, the lowest eigenvalue of G<g is greater than that 



of ju(5 r ) 1 Gs r - This inequality then yields (27), since these eigenvalues are 



IGv^ 1 ^ 1 and //(S' r .)||G s 1 || 2 1 , respectively. □ 



We now need to compute the Gram matrix for the canonical basis of Hl . 
This basis is described in |24l Chapter IX, §3.6] and consists of spherical 
harmonics. Let £, k\, . . . , kd-i be integers satisfying I > k\ > k% > ■ ■ ■ > 
kd-i > 0, and take K := (ki, . . . , ±/^_i). A spherical harmonic of degree 
£ O p. 466]) will be denoted by Yg(9i, ... ,9d)- The angles are the usual 
ones from spherical coordinates in IR rf+1 (cf. |24l p. 435]). The basis for LT^ 
is then the set of all Y^, < £ < L. The entries in the Gram matrix are 
[Gs r ](l,K),(t',K') = ( Y K' Y K')s r - Following the argument in [Ml Chapter IX, 
§3.6], one may show that 



d-l 



= B 1jK Bv jK 5 k ,k> C t \ '(cos e)C e 2 _ k \ L (cos0)sin 2Kl+a - l 9de, (28) 
J o 

where C^(t) is the Gegenbauer polynomial of degree n and type s, and Be fa 
is a normalization factor. In the case where d = 2 and L = 1, Gs T is 4 x 4 
and has six non-zero entries, 



£(0,0), (o,o) — |(1 



cos r 



G ! (o,o),(i,o) = G(i j o),(o,o) = XC 1 -cosr)(l + cosr), 

G (l,0),(l,0) = |(l-cosr)(l + cosr + cos 2 r), 
G r (i,±l),(i,±l) = i(l -cosr) 2 (2 + cosr). 

Since fJ,(S r ) = 2-/r(l — cosr), the formulas for the entries above imply that 
Gs r l n(S r ) is a polynomial in cosr. In fact, a straightforward calculation 
shows t hat the minimum eigenvalue of this matrix is r 4 /(256-7r) + 0(r e ). 
Lemma Eg] then implies that \\G^ X \\ 2 < n(S r )\\G s x \\ 2 = 2567rr" 4 + 0(r- 2 ). 
A less precise, but similar result, holds in the general case. 



Lemma 6.4. Under the assumptions of Lemma 6.3, for general L > 0, 
d > 2 and r sufficiently small, there is an integer 1 = l(L, d) > L and a 
constant C = C(L, d) > such that \\G7 {[2 < Cr~ 2L . For L = 1 and d = 2, 
we may take 1 = 2. 
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Proof. From the expression in ( 28 ) for the entries in Gs r , we see that each 
of them is entire in r and has a zero of order d or greater at r = 0. In 
addition, fj,(S r ) is also entire in r and has a zero of order d. It follows that 
the matrix G{r) = Gs r / n(S r ) is entire, even in r. and for real r, it is real, 
self adjoint and positive semi definite. (In fact, for d even, it is a polynomial 
in cosr.) In addition, the 2x2 block in G(0) corresponding to k\ = 0, 
t = 0, 1, is rank 1 and therefore has as an eigenvalue; consequently, G(0) 
also has as an eigenvalue - it's lowest, in fact. As Rellich [191 Pg- 91] 
shows, the eigenvalues of G(r) are analytic functions of r. For r > 0, these 
eigenvalues are proportional to those of the Gram matrix Gs r and therefore 
must be positive. None of these eigenvalues are identically 0. In particular, 
the eigenvalues splitting off from the eigenvalue of G(0) are not identically 
0. As functions of r they thus have a zero of finite order at r = 0; the order 
is an even integer because G{r) is even in r. The smallest eigenvalue then 
behaves like A m j n (r) = r 2t (ao + C(r 2 )), where ao > 0, i > is an integer, 
and r is sufficiently small. Furthermore, from (28) we see that the diagonal 
entry, with I = £' = k± = L, is 0(r 2L ). since this bounds the minimum 
eigenvalue from abo ve, we must have 2t > 2L, so t > L. The result then 
follows from Lemma |6\3| and the observation that M-S'rOIIG^ 1 !^ = A~* n . The 
calculation for L = 1 and d = 2 was done above. □ 



6.3 Local Lagrange Bases 

We now turn to the local Lagrange basis. Recall that the function G 
S m (E), with the kernel representation 

m 2 

Xi = J2 A U k m(;0 + J2' b ^J G S m (~), (29) 

is a local Lagrange function centered at £ if it satisfies x?|t = e £> where 
e^(() = 5^^; that is, is the vector (1, 0, . . . , 0) T . Since \£, £ S m (E), 
the vector = ,f)feT is in the constraint space. This vector and the 
coefficients bj then satisfy x?|t = e £ ; Of course, the (full) Lagrange function 

X£ = 2~^ e H ^£,C K ('> + Y?jLibj<f>j restricted to T also satisfies x^lr = 
e^. Consequently, 25 £ := — X( satisfies 25^|x = 0. We can rewrite this 
difference as 25^ = X£ — + X£ ~ X£ (with xg the truncated basis function 
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introduced in the last section). It follows that 



#J 



Pi 



XC-X5 



Evaluating this on T then gives the system Kja + Qft + (xg — X?)|t = 0. 
By linearity, it is clear that <I>*Q! = 0. Finally, letting y = (xs, ~ Xc)|x> we 
arrive at the system, 



K T 3>\ fa 
0/1/3 



(30) 



Proposition 5.2 applies to (30), with H replaced by T; thus, noting that 
||y||^(T) < \\X£ - XtWoo, and writing J = (1, . . . ,m 2 ), we see that 

IHk(r) < $~ 1 \/Wt\\xs, - xelloo 

< 2||A: m || 00 ||G T 1 || 1 /^- 1 (#T) 3 / 2 ||x C -X f lloo 
From this we obtain these inequalities: 

Wxt ~ Xglloo < \\km\\oo\/W^\\oi\U 2 rc) + mCm\\f3\\e 2 (j), C m = max ||^||oo 

j<m' ! 

< ||A: m ||oo#T^- 1 + 2mC m ^#T||G T 1 ||^ - 
Moreover, using ||x£ - X£lU < ~ + ||X£ - Xflloo, we see that 



X^-xdU <2||A: m || 00 #T ? ?- 1 ( l + 2mC m d#T\\G^\\ ) ||x£ - Xclloo- (31) 



-f—i i 



Finally, from Proposition 6.1, if K > (2m + 2(j,)/v, it is easy to see that this 
holds: 

llfc - xdU < C^|^ (1 + mh-») h K »- 2m -» < C ^^-h Kv - 2m - 2 » 



r) 



i) 



(32) 

To proceed further, we need to estimate Such estimates are known for 
surface splines in the Euclidean case [181 §6]- Simply repeating the proofs of 
|18} Corollary 2.2] and [18} Theorem 2.4] for a set of points in M 3 restricted 
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to S 2 yields the desired estimate. For the collocation matrix associated with 
k m and S, we have 

> Cq 2m ~ 2 , (33) 

where C depends only on mj^] Thus, for k m , we have ||x^ — X^Woo < 
(Jh Ku ~ Am+2 ~ 2fl , where the constant X has to be increased slightly to absorb 
| log /i| 2 . With this in mind we have the following result, whose proof, being 



similar to Proposition 6.1, we omit 



Theorem 6.5. Let the notation and assumptions of Theorem 5.3 hold. Sup- 
pose that K > is chosen so that K > 4m ~2+ 2 M and, for each £ G S, T(£) : = 
EnB(^, Kh\ logh\). Ifx$ * s a local Lagrange function for T(£) centered at£, 
then set £eH is a basis for S m {z>) . Moreover, with J := X^— 4m+2 — 2/x, 

||x ? -X 5 lloo<C/i J , (34) 
|fc(x)| < C(l + dist(x,e)// l )- J . (35) 

Furthermore, when J > d, the set {x^} is L p stable: there are Ci,C% > 
for which 

Cig 2/p ||a||, p(H) < ||E ?eS a?xdL p(s2) < <V /p l|a||*<E). (36) 
Quasi-interpolation: It follows that the operator 

provides L^ convergence at the same asymptotic rate as interpolation /=. 
Indeed, 

\I~f(x)-Q s f(x)\<J2\Mx)-X&)\\m)\ < cq-^fuh^-^- 2 ^ 

< C\\fUh 2m 



4 For any d-dimensional sphere or projective space and any conditionally positive def- 
inite polyharmonic kernel with associated polynomial operator C m = Q(— A), where 
Q is a polynomial of degree m, the coefficients in the expansion for k m are given by 
k(j) — Q(Xj)^ 1 , j £ J . For large \j, all of these have the asymptotic behavior A~ m , 
which is the same as that of the coeffici ents for the m-d thin-plate spline. This implies 



that the matrix P ± KrP ± in Proposition 5.2 (here, 3 — > T) will have a lowest eigenvalue 



value that is, up to a constant multiple, dependent only on m and d. Consequently, the 
bound i? > Cq 2m ~ d holds for all k m associated with C m in dimension d. 
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provided that K > 6m +^+ 2 1 n j s shown in Corollary 5.9] that restricted 
surface spline interpolation exhibits ||ia/ — /||oo — Ch u for / G C 2m (E> 2 ) 
when a = 2m and for / E B^ ^(S 2 ) for a < 2m. So Qs has the same rate 
of approximation (without needing to solve a large system of equations) . 

Constructing basis functions in terms of N. Given a set of scat- 
tered points, it may be desirable to use N as the basic parameter instead 
of h. Therefore we wish to express the number of nearest neighbors needed 
as a function of the total cardinality N instead of those within a Kh log h 
neighborhood. Considering a cap B(a,r), a simple volume argument gives 

36 /r x 2 



#(%r)nS)<-M . (37) 

Indeed, one arrives at this bound by first considering caps of radius q around 
each node in B(a,r). If q is small enough, say q < r, then at least 1/3 the 
volume of each cap will be contained in B(a,r). Using this and a Taylor 
expansion of the volume formula 2tt(\ — cos(q)) leads to ([37]). Thus, the 
greatest number of points in a cap of radius Kh\ log h\ is |yp 2 (K log(l//i)) 2 . 
Also, it is not hard to show that 2h~ 2 < N, and hence it follows that the 
number of points is bounded by ffp 2 ( : f log(iV)) 2 = Ti(p^) 2 (l°g(-^0) 2 > an d 
it suffices to take for T the nearest jj(pK) 2 (log(iV)) 2 neighbors. 

The constants v and K. Before we turn to a discussion of precondi- 
tioning, we wish to comment on the constants v and K above. These two 
constants come into play in a crucial way in many of our estimates. 



The decay constant v first comes up in the proof of Theorem 5.3 (Al- 
though we do not mention it in the theorem, the proof produces two different 
decay constants: vl and vc, the former for the Lagrange function and the 
latter for the coefficients.) The estimate for v is, of course, a lower bound 
on the decay constant itself; it is independent of p, but weakly dependent 
on m. Because of the nature of such estimates, it is very likely that they are 
much lower than f ac tuai- How z^ actua i behaves as a function of p is an open 
question. 

There is another open question concerning K. We know that it must be 
bounded below by 4m ~ 2+2/J . Thus a better estimate on v would produce 
a better lower bound on K. This in turn means using smaller caps and 
fewer points in constructing the local Lagrange interpolant - i.e., giving it 
a smaller "footprint." On the other hand, the larger we make K the better 
the approximation to xs, we § e t- Since K can be made as large as we please, 
the question then becomes this: What is an optimal choice for Kl Indeed, 
what does the term optimal mean here? 
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7 Preconditioning with local Lagrange functions 

In this section we illustrate how the local Lagrange functions can also be 
used as an effective preconditioner for linear systems associated with in- 
terpolation using the standard restricted spline basis. Our focus is on the 
restricted surface spline &2 (i.e. the restricted thin plate spline), for which 
the interpolant to /L in the standard basis takes the form 



Ses j=i 



(38) 



where 4>j are a basis for the spherical harmonics of degree < 1. We note 
that this interpolant can also be written with respect to the local Lagrange 
basis for ^(H) as 



(39) 



see Section [3] for the details on constructing this basis. 

Using the properties of the local Lagrange basis, we can write the linear 



system for determining the interpolation coefficients ag in (39) as: 



[Kg $] 



(40) 



where (K H )y = i,j = 1,...,N, and = i = l,...,N, 

j = 1, . . . , 4. The matrix Ax is a N-by-N sparse matrix where each column 
contains n = M(\ogN) 2 entries corresponding to the values of the interpo- 
lation coefficients for the local Lagrange basis in Q. The matrix Cx is 
a 4-by-A^ matrix with each column containing the values of the interpolation 
coefficients cgj in ([6]). With the linear system written in this way, one can 
view the matrix [Ax Cy] T as a right preconditioner for the standard kernel 



interpolation matrix. Once a is determined from ( 40 ) , we can then find the 



interpolation coefficients and Cj in (38) from 



[a c] T = [A T C T ] T a. 



(41) 



If the local Lagrange basis decays sufficiently fast then the linear system 
(40) should be "numerically nice" in the sense that the matrix K=Ax + 



<I>Cx should have decaying elements from its diagonal and should be well 
conditioned. As discussed in the previous section, the decay is controlled by 
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Figure 5: Interpolated relative vorticity from a numerical simulation of the 
shallow water wave equations on the N = 163842 icosahedral node sets. 
The original values for the relative vorticity come from [H] and have been 
interpolated to a regular 300 x 600 latitude-longitude based grid using the 
restricted kernel spline k2(x, a) = (1 — x ■ a) log(l — x ■ a). The interpolation 



coefficients were computed using GMRES on the preconditioned system ( 40 ) . 



the number of nearest neighbors n used in constructing each local Lagrange 
function and that n = M(log iV) 2 . In the experiments below, we found that 
choosing n = 7[(log iV) 2 /(log 10) 2 )] = 7[(log 10 N) 2 ] gave very good results 
over several decades of N. 



To solve the preconditioned linear system ( 40 ) we will use the generalized 



minimum residual method (GMRES) [21J. This is a Krylov subspace method 
which is applicable to non-symmetric linear systems and only requires com- 
puting matrix-vector products. Each matrix-vector product involving the 
preconditioner matrix [Aj- Cx] T requires 0(iV(log -/V) 2 ) operations, while 
each matrix-vector product involving [K= 3>] requires 0(N 2 ) operations. 
However, Keiner et al. have shown that this latter product can be done 
in 0{N log N) operations using fast algorithms for spherical Fourier trans- 
forms |13j . As we are primarily interested in the exploring the effectiveness 
of the local Lagrange basis as a preconditioner, we have not used these fast 
algorithms in the results below. In a follow up study, we will investigate 
these fast algorithms in combination with the preconditioner in much more 
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detail. 

For the first numerical tests we use icosahedral node sets 3 C S 2 of in- 
creasing cardinality. These were chosen because of their popularity in com- 
putational geosciences (see, for example, [9] [22], [20] [15] ) where interpolation 
between node sets is often required. The values of / were chosen to take on 
random values from a uniform distribution between [—1,1]. Table [2] displays 
the number of GMRES iterations to compute an approximate solution to 



the resulting linear systems (40) for various N and different tolerances. As 
we can see, the number of iterations is small and stays relatively constant 
as ./V increases. 







Number GMRES iterations 


N 


n 


tol = 10~ 6 


tol = 10~ 8 


2562 


84 


7 


5 


10242 


119 


5 


7 


23042 


140 


6 


7 


40962 


154 


5 


7 


92162 


175 


6 


8 


163842 


196 


5 


7 



Table 2: Number of GMRES iterations required for computing an approx- 



imate solution to (40) using icosahedral node sets of cardinality N. Here n 
corresponds to the number of nodes used to construct the local basis and 
tol refers to the tolerance on the relative residual in the GMRES method. 
The right hand side was set to random values uniformly distributed between 
[—1,1] and the initial guess for GMRES was set equal to the function values. 



For the final numerical experiment, we use the above technique to in- 
terpolate a field taken from a numerical simulation on the icosahedral node 
sets to a regular latitude-longitude grid. As mentioned above, this is often 
necessary for purposes of comparing solutions from different computational 
models, plotting solutions, or coupling different models together. The data 
we use comes from [8] and represents the relative vorticity of a fluid described 
by the shallow water wave equations on the surface of a rotating sphere. The 
initial conditions for the model lead to the development of a highly nonlin- 
ear wave with rapid energy transfer from large to small scales, resulting in 
complex vortical dynamics. The numerical solution was computed on the 
N = 163842 node set and we interpolated it to a regular 300 x 600 latitude- 
longitude based grid. Figure [5] displays the resulting interpolated relative 
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vorticity from the simulation at time t = 6 days. The figure clearly shows 
that the complex flow structure has been maintained after the interpola- 
tion. As in the numerical examples above, the approximate solution to (40) 
with this data was obtained in 7 iterations of the GMRES method using a 
tolerance of 10~ 8 . 
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